Metode numerice - Aplicatii

Lucrarea 5.   Rezolvarea sistemelor de ecuatii liniare prin factorizare Crout si Choleski

Tema A - Factorizarea Crout                Tema B - Factorizarea Cholesky

        Prin definitie, descompunerea unei matrice A(n,n) intr-un produs de doua matrice:

A = L * R
unde L(n,n) este o matrice inferior triunghiulara, iar R(n,n) este o matrice superior triunghiulara, se numeste factorizare LR a matricei A. De exemplu, pentru o matrice A de dimensiuni 4 * 4  factorizarea LR are forma:
        In cazul folosirii factorizarii, rezolvarea unui sistem de ecuatii liniare A * x = b se reduce la rezolvarea succesiva a doua sisteme triunghiulare:
A * x = b    <==>    ( L * R ) * x = b    <==>    L * ( R * x ) = b
Mai intai se rezolva sistemul inferior triunghiular:
L * y = b
in raport cu y (prin substitutie inainte) si apoi sistemul superior triunghiular:
R * x = y
in raport cu x (prin substitutie inapoi) , necunoscutele principale ale problemei.
        Aplicarea procedurilor compacte de factorizare (pe loc, in matricea A) presupune determinarea a n^2+n elemente ale matricelor L si R, folosind cele n^2 ecuatii care se pot scrie. Pentru ca sistemul de ecuatii astfel obtinut sa admita solutie unica este necesara precizarea a priori a n din cele n^2+n  necunoscute. Metodele de factorizare utilizate aleg ca necunoscute dependente, ale caror valori se precizeaza a priori, cele n elemente diagonale ale uneia din matricele L sau R. Se disting astfel doua tipuri de factorizari:

    (i) factorizarea Crout, care impune matricea R cu diagonala unitate;
    (ii) factorizarea Doolitle, care impune matricea L cu diagonala unitate.

        Daca matricea A este simetrica si pozitiv definita, cele doua matrice de factorizare satisfac relatia, astfel incat factorizarea devine:

cunoscuta sub numele de factorizare Cholesky.

Tema A - Factorizarea Crout

        Impunand elementele diagonale ale matricei R egale cu unitatea, factorizarea Crout determina restul elementelor necunoscute din matricele L si R prin  rezolvarea unui sistem format din n^2 ecuatii, cu tot atatea necunoscute. Ecuatiile care formeaza acest sistem se deduc din relatia de definitie A = L * R si forma lor depinde de raportul in care se afla indicii i si  j ai unui element a_ ij din matricea A. Putem avea urmatoarele trei cazuri:

Cele trei ecuatii de mai sus pot fi scrise compact sub forma:
           (**)
Tehnica factorizarii Crout nu realizeaza rezolvarea propriu-zisa a acestui sistem de ecuatii, ci determina necunoscutele _ij  si _ij printr-o aranjare convenabila a ecuatiilor.

        Astfel, vom considera ca au fost determinate elementele nenule de pe primele p-1 coloane din matricea L si primele p-1 linii din matricea R si vom particulariza relatia (**) pentru indicii i si  j avand succesiv valoarea p.
        Pentru determinarea elementelor nenule de pe coloana p a matricei L, in relatia (**) se impune j=p si se individualizeaza elementele acestei coloane. Deoarece pe coloana p a matricei L singurele elemente nenule_ip corespund unor indici i=p,...,n, limita de sumare va fi egala cu p, astfel incat se poate scrie:

Am admis insa ca toate elementele de pe primele p-1 coloane din matricea L (_im; m=1,...,p-1), respectiv de pe primele p-1 linii din matricea R (_mp; m=1,...,p-1) au fost determinate in prealabil. De asemenea, _pp = 1 (ipoteza initiala a metodei). In aceste conditii, relatia de mai sus permite calculul tuturor elementelor nenule de pe coloana p a matricei L:
           ( I )
        Pentru determinarea elementelor nenule de pe linia p a matricei R, in relatia (**) se impune i=p si se individualizeaza elementele acestei linii. Totodata, deoarece pe linia p a matricei R elementele nenule necunoscute _pj corespund unui indice  j=p+1,...,n (deoarece _pp=1 este cunoscut), limita de sumare este egala cu p, astfel incat relatia (**) capata forma:
Si in acest caz toate elementele de pe primele p-1 coloane din matricea L (_pm; m=1,...,p-1), respectiv de pe primele p-1 linii din matricea R (_mj; m=1,...,p-1) au fost determinate in prealabil. De asemenea, elementul diagonal _pp din matricea L a fost determinat anterior. Astfel, relatia de mai sus permite calculul elementelor nenule de pe linia p a matricei R:
       ( II )
Aplicarea relatiilor ( I ) si ( II ) pentru toate valorile lui p=1,...,n permit determinarea tuturor elementelor nenule din matricele L si R, realizand astfel factorizarea Crout a matricei A.

Aceste relatii evidentiaza faptul ca la calculul elementelor  nenule de coloana p a matricei L, respectiv de pe linia p a matricei R, se folosesc numai elemente de pe coloana, respectiv linia p a matricei A. Ca urmare, noile valori calculate _ip  si _pj  pot fi suprapuse peste elementele a_ip, respectiv a_pj, inlocuindu-le pe acestea din urma. Astfel, factorizarea Crout se poate desfasura pe loc, in matricea A, ne fiind necesara ocuparea de memorie auxiliara pentru matricele L si R.


Algoritmul 1 - Rezolvarea sistemelor de ecuatii liniare cu factorizarea Crout
  1. Definirea sistemului de ecuatii: rangul n, matricea coeficientilor A, vectorul termenilor liberi b.
  2. Factorizarea LR a matricei A dupa metoda Crout :
    1. Initializarea numarului coloanei / liniei din matricele L si R pentru care se calculeaza elementele nenule: p <-- 1;
    2. Determinarea elementelor subdiagonale, inclusiv a elementului diagonal de pe coloana p a matricei L (relatia I):
    3. Pivotarea partiala pe coloana p;
    4. Determinarea elementelor situate la dreapta diagonalei principale , pe linia p a matricei R (relatia II ):
    5. Se trece la tratarea urmatoarei coloane/linii din matricele L si R: p <--  p+1. Daca p <= n se revine la pasul 2.2. Daca p > n se trece la pasul 3.
  3. Rezolvarea sistemului inferior triunghiular L * y = b, prin substitutie inainte, folosind elementele diagonale si subdiagonale din matricea A si memorarea solutiei y in vectorul b.
  4. Rezolvarea sistemului superior triunghiular R * x = y = b prin substitutie inapoi, folosind elementele supradiagonale din matricea A si memorarea solutiei x in vectorul  b.

Pentru modul in care factorizarea Crout poate fi folosita la inversarea unei matrice, se poate consulta cartea "Calcul numeric cu aplicatii in Turbo Pascal".
 

Tema B - Factorizarea Cholesky

        Daca o matrice patrata A este simetrica si pozitiv definita, ei i se poate aplica o factorizare speciala, mai eficienta din punct de vedere numeric, numita factorizare Cholesky. Simetria inseamna satisfacerea egalitati  (T indica transpunerea) sau a_ ij = a_ ji pentru toti indicii i,j=1,...,n. O matrice A se numeste pozitiv definita daca inegalitatea:

este satisfacuta pentru orice vector x, iar anularea produsului matriceal are loc numai atunci cand x este identic cu vectorul nul. Proprietatea definirii pozitive a matricei A joaca un rol esential in stabilitatea numerica a calculelor efectuate in cursul factorizarii Cholesky.
        Daca cele doua proprietati sunt satisfacute, factorizarea Cholesky descompune matricea A intr-un produs de doua matrice triunghiulare de tipul LR, cu proprietatea remarcabila  :
Pentru acest tip de factorizare, se pot scrie (n^2+n)/2 ecuatii independente care contin tot atatea necunoscute (elementele nenule din matricea L). In consecinta, factorizarea Cholesky a unei matrice este unica, nefiind necesara precizarea a priori a valorilor unor necunoscute.
Pentru stabilirea relatiilor generale de calcul dupa care se desfasoara factorizarea Cholesky, expresia se expliciteaza in functie de elementele matricelor A si L:
        (***)
        In ipoteza determinarii prealabile a elementelor de pe primele  j-1 coloane ale matricei L, relatia (***) se expliciteaza pentru a evidentia elementul diagonal  :
si se face particularizarea i = j:
Deoarece toti termenii _ jm ( m=1,...,j-1 ) au fost calculati in prealabil, se poate determina elementul diagonal  :
        ( III )
dupa care, se calculeaza restul elementelor de pe coloana j a matricei L:
        ( IV )
Relatiile ( III) si ( IV ) stau la baza factorizarii matricei A dupa metoda Cholesky. Aplicarea lor repetata, pentru toti indicii  j=1,...,n, permite determinarea celor n coloane ale matricei L.

        Si in cazul factorizarii Cholesky exista posibilitatea aplicarii metodei descrise in forma compacta, pe loc in matricea A. Pe de alta parte,  este posibila memorarea intr-o singura matrice atat a matricei originare A, cat si a matricei de factorizare L. Astfel, forma finala a matricei A are urmatoarea structura: pe diagonala si deasupra acesteia se vor gasi elementele originare ale matricei A; sub diagonala se afla elementele matricei L, iar diagonala acesteia este memorata separat intr-un vector.
        O particularitate a factorizarii Cholesky se refera la aplicarea tehnicilor de pivotare. In acest sens, daca matricea A satisface conditiile pentru a admite o factorizare Cholesky (este simetrica si pozitiv definita), ea va fi întotdeauna nesingulara. Prin urmare riscul impartirii la un element nul in relatia ( IV ) nu exista si nici erorile de rotunjire nu pot afecta semnificativ rezultatele calculelor. Prin urmare, nu mai este necesara aplicarea pivotarii. Mai mult decat atat, aplicarea tehnicilor de pivotare nici nu este permisa deoarece:

  1. pivotarea partiala strica simetria matricei;
  2. pivotarea completa pastreaza simetria, dar conduce, de regula, la pierderea proprietatii de "definire pozitiva" a matricei A.
In ambele cazuri se incalca conditiile care asigura existenta factorizarii Cholesky.


Algoritmul 2 - Rezolvarea sistemelor de ecuatii liniare cu factorizarea Cholesky
  1. Definirea sistemului de ecuatii: rangul n, matricea coeficientilor A, vectorul termenilor liberi b.
  2. Verificarea simetriei matricei A. Daca matricea nu este simetrica se lanseaza un mesaj de eroare si se intrerupe algoritmul;
  3. Factorizarea Cholesky pe loc in matricea A:
    1. Initializarea coloanei curente j din matricea L pentru care se determina elementele nenule: j <-- 1;
    3. Calculul expresiei de sub radical pentru determinarea elementului diagonal_ jj :
      si verificarea semnului acesteia. Daca S < 0 se transmite un mesaj de eroare si se intrerupe algoritmul. Daca S >= 0 se trece la pasul 3.iii;
    4. Determinarea elementului diagonal _ jj si memorarea lui intr-un vector distinct:
      D_ j <--;
    5. Determinarea elementelor subdiagonale de pe coloana  j a matricei L si memorarea lor in matricea A:
    6. Se trece la tratarea urmatoarei coloane din matricea L : j <-- j+1. Daca  j <= n se revine la pasul 3.2. Daca j > n se trece la pasul 4;
  4. Rezolvarea sistemului inferior triunghiular L * y = b prin substitutie inainte si memorarea solutiei in vectorul b;
  5. Rezolvarea sistemului superior triunghiular * x = y = b prin substitutie inapoi si memorarea solutiei in vectorul b.

Pentru detalii suplimenatare privind factorizarea Cholesky, se poate consulta cartea "Calcul numeric cu aplicatii in Turbo Pascal".
 

   Aplicatii - Lista lucrarilor de laborator